#include "header.h"

/* --------------------------------------------------------------------- */
void geometric(struct node* p){
/* --------------------------------------------------------------------- */
int i,j,k,l,neigh,LSIZE,iMaxCond,iMaxCondWeighted,status;
double dx,dy,dz,dx2,dy2,dz2,dr2,dxdy,det,ds,ds2,weight[50],cond,DIST;
double wdx2,wdy2,wdr2,wdxdy,wdet,wcond,wt;
double sdx,sdy,**Lmat,Amat[50][6],b[6];

printf("Finding geometric coefficients...");
MaxCond = 1.0;
MaxCondWeighted  = 1.0;
iMaxCond = 0;
iMaxCondWeighted = 0;
for(i=0;i<no_grid_points;i++){
		dx2 = dy2 = dz2 = 0.0;
            for(j=0;j<p[i].total_conn;j++){
                neigh = p[i].Conn[j];
                dx = p[neigh].x - p[i].x;
                dy = p[neigh].y - p[i].y;
                dz = p[neigh].z - p[i].z;

                ds2   = (dx*dx + dy*dy + dz*dz);
                wt    = WEIGHT(ds2, 1.0);

		    dx2 += dx*dx;
		    dy2 += dy*dy;
		    dz2 += dz*dz;
        }
		dr2      = dx2 + dy2 + dz2;
		p[i].ds2 = dr2/p[i].total_conn;
		if(p[i].ds2 <= 0.0){
			printf("Distance is negative for node %d.\n",i);
			exit(0);
			}
 
} 

/* Find minimum length for each connectivity which used in
   time-step calculation */
MinEdge = 1000000.0;
MaxEdge = 0.0;
DsVenkat = 0.0;
for(i=0;i<no_grid_points;i++){
	DIST = 100000.0;
	DsVenkat += p[i].ds2;
	for(j=0;j<p[i].total_conn;j++){
		neigh = p[i].Conn[j];
		dx    = p[neigh].x - p[i].x;
		dy    = p[neigh].y - p[i].y;
		dz    = p[neigh].z - p[i].z;
		ds    = sqrt(dx*dx + dy*dy + dz*dz);
		if(ds == 0.0){
			printf("Distance zero between node  %d  and  %d\n",i,neigh);
			exit(0);
			}
		DIST = min(DIST, ds);
		}
	p[i].dsm = DIST;
	MinEdge  = min( MinEdge, DIST );
	MaxEdge  = max( MaxEdge, DIST );
	}

DsVenkat = sqrt( DsVenkat/no_grid_points );

/* Find coefficients for derivative */
LSIZE = 3;
Lmat = (double**)calloc(LSIZE,sizeof(double*));
for(i=0;i<LSIZE;i++) Lmat[i] = (double*)calloc(LSIZE,sizeof(double));
for(i=0;i<no_grid_points;i++){
	b[0] = b[1] = b[2] = b[3] = b[4] = b[5] = 1.0;
	for(j=0;j<p[i].total_conn;j++){
		neigh = p[i].Conn[j];
		dx    = p[neigh].x - p[i].x;
		dy    = p[neigh].y - p[i].y;
		dz    = p[neigh].z - p[i].z;
		ds2 = dx*dx + dy*dy + dz*dz;
		weight[j] = WEIGHT(ds2,p[i].ds2);
		Amat[j][0] = dx;
		Amat[j][1] = dy;
		Amat[j][2] = dz;
		}
	for(j=0;j<LSIZE;j++)
		for(k=0;k<LSIZE;k++){
			Lmat[j][k] = 0.0;
			for(l=0;l<p[i].total_conn;l++)
				Lmat[j][k] += Amat[l][j] * weight[l] * Amat[l][k];
			}
	status=Gaussj(Lmat,LSIZE,b);
	if(status==0){
		printf("Fatal: could not determined coefficients\n");
		printf("       Node = %d\n", i);
		exit(0);
		}
	for(j=0;j<p[i].total_conn;j++){
		p[i].X[j] = p[i].Y[j] = p[i].Z[j] = 0.0;
		for(k=0;k<LSIZE;k++){
			p[i].X[j] += Lmat[0][k] * weight[j] * Amat[j][k];
			p[i].Y[j] += Lmat[1][k] * weight[j] * Amat[j][k];
			p[i].Z[j] += Lmat[2][k] * weight[j] * Amat[j][k];
			}
		}
	}

for(i=0;i<LSIZE;i++) free(Lmat[i]); free(Lmat);
printf("Done\n");
}
	
